#ifndef _BLES7_H_
#define _BLES7_H_

#include "Equations/TwoPhaseEquation.hpp"

namespace Tuna {

  template <typename T, int Dim> class BLES7;

  template <typename T>
  struct Typeinfo <BLES7<T, 1> > {
    typedef T prec_t;
    enum { Dim = 1 };
  };

  template <typename T>
  struct Typeinfo <BLES7<T, 2> > {
    typedef T prec_t;
    enum { Dim = 2 };
  };

  template <typename T>
  struct Typeinfo <BLES7<T, 3> > {
    typedef T prec_t;
    enum { Dim = 3 };
  };
  
  template<typename Tprec, int Dim>
  class BLES7 : public TwoPhaseEquation<BLES7<Tprec, Dim> >
  {

    typedef TwoPhaseEquation<BLES7<Tprec, Dim> > TP_BLES7;

    using GeneralEquation< TP_BLES7 >::aE;
    using GeneralEquation< TP_BLES7 >::aW;
    using GeneralEquation< TP_BLES7 >::aN;
    using GeneralEquation< TP_BLES7 >::aS;
    using GeneralEquation< TP_BLES7 >::aF;
    using GeneralEquation< TP_BLES7 >::aB;
    using GeneralEquation< TP_BLES7 >::aP;
    using GeneralEquation< TP_BLES7 >::sp;
    using GeneralEquation< TP_BLES7 >::dx;
    using GeneralEquation< TP_BLES7 >::dy;
    using GeneralEquation< TP_BLES7 >::dz;
    using GeneralEquation< TP_BLES7 >::dt;
    using GeneralEquation< TP_BLES7 >::bi;
    using GeneralEquation< TP_BLES7 >::ei;
    using GeneralEquation< TP_BLES7 >::bj;
    using GeneralEquation< TP_BLES7 >::ej;
    using GeneralEquation< TP_BLES7 >::bk;
    using GeneralEquation< TP_BLES7 >::ek;
    using GeneralEquation< TP_BLES7 >::applyBoundaryConditions1D;
    using GeneralEquation< TP_BLES7 >::applyBoundaryConditions2D;
    using GeneralEquation< TP_BLES7 >::applyBoundaryConditions3D;

    using TP_BLES7::S;
    using TP_BLES7::p;
    using TP_BLES7::phi_0;
    using TP_BLES7::Srw;
    using TP_BLES7::Sro;
    using TP_BLES7::mu_w;
    using TP_BLES7::mu_o;
    using TP_BLES7::k;
    using TP_BLES7::porosity;

  public:
    typedef Tprec prec_t;
    typedef typename TunaArray<prec_t, Dim >::huge ScalarField;
    
    BLES7() : TwoPhaseEquation<BLES7<prec_t, Dim > >() { }    
    ~BLES7() { };
    
    inline bool calcCoefficients1D(); 
    inline bool calcCoefficients2D();
    inline bool calcCoefficients3D();
    inline void printInfo() { std::cout << " BLES7 "; }
  };

/*
 *  Lineal for realtive permeability , QUICK for Sw
 */
template<typename Tprec, int Dim>
inline bool BLES7<Tprec, Dim>::calcCoefficients1D () 
{
    static prec_t Sw_e, Sw_w;

    // Lineal 
    static prec_t mult = k * dt / 
      ( porosity * dx * dx * (1 - Srw - Sro) * mu_w );
      
    aE = 0.0; aW = 0.0; aP = 0.0; sp = 0.0;

    for (int i = bi; i <= ei; ++i) {      

      // UpwindQ
      if ( p(i+1) >= p(i) )
	if      ( i == ei   ) Sw_e = phi_0(i+1);
	else if ( i == ei-1 ) Sw_e = phi_0(i+1);
	else Sw_e = phi_0(i+1) + 0.125 * ( 3 * phi_0(i) - 2 * phi_0(i+1) - phi_0(i+2) );
      else
	if ( i == bi ) Sw_e = phi_0(i);
	else Sw_e = phi_0(i) + 0.125 * (3 * phi_0(i+1) - 2 * phi_0(i) - phi_0(i-1) );
      
      if ( p(i-1) >= p(i) ) 
	if      ( i == bi   ) Sw_w = phi_0(i-1);
	else if ( i == bi+1 ) Sw_w = phi_0(i-1);
	else Sw_w = phi_0(i-1) + 0.125 * (3 * phi_0(i) - 2 * phi_0(i-1) - phi_0(i-2) );
      else
	if ( i == ei ) Sw_w = phi_0(i);
	else Sw_w = phi_0(i) + 0.125 * (3 * phi_0(i-1) - 2 * phi_0(i) - phi_0(i+1) );

      // Lineal
      aE (i) = (Sw_e - Srw) * mult;
      aW (i) = (Sw_w - Srw) * mult;
      aP (i) = aE (i) + aW (i);      
    }

/* ----- Boundary conditions ----- */

    // I'm using an explicit scheme, so the boundary conditions are
    // different than in the implicit case.

    // Dirichlet right side
    aP(ei) += aE(ei);
    sp(ei) = 2 * aE(ei) * p(ei+1);
    aE(ei) = 0;

    // Neumann left side
    // Remember that Sw_e = 0.8 on the boundary, that is
    // the reason why sp(bi) has the form below.
    aP(bi) -= aW(bi) ;
    sp(bi) = aW(bi) * dx * ( 3.4722e-7 * mu_w / k) ;
	//   sp(bi) = aW(bi) * dx * ( 3.00e-4 * mu_w / k) ;
    aW(bi) = 0;
    
    return 0;
}


/*
 *  Lineal for realtive permeability , Upwind for Sw
 */
template<typename Tprec, int Dim>
inline bool BLES7<Tprec, Dim>::calcCoefficients2D () 
{
  static prec_t Sw_e, Sw_w, Sw_n, Sw_s;

    // Lineal 
  static prec_t multx = k * dt / 
    ( porosity * dx * dx * (1 - Srw - Sro) * mu_w );

  static prec_t multy = k * dt / 
    ( porosity * dy * dy * (1 - Srw - Sro) * mu_w );

  aE = 0.0; aW = 0.0; aN = 0.0; aS = 0.0; aP = 0.0; sp = 0.0;
  
  for (int i =  bi; i <= ei; ++i)
    for (int j = bj; j <= ej; ++j) {  

      // Upwind
      if ( phi_0(i+1, j) >= phi_0(i, j) ) Sw_e = phi_0(i+1, j);
      else                                Sw_e = phi_0(i  , j);
      if ( phi_0(i-1, j) >= phi_0(i, j) ) Sw_w = phi_0(i-1, j);
      else                                Sw_w = phi_0(i  , j);

      if ( phi_0(i, j+1) >= phi_0(i, j) ) Sw_n = phi_0(i, j+1);
      else                                Sw_n = phi_0(i, j);
      if ( phi_0(i, j-1) >= phi_0(i, j) ) Sw_s = phi_0(i, j-1);
      else                                Sw_s = phi_0(i, j);

      // Lineal
      aE (i, j) = (Sw_e - Srw) * multx;
      aW (i, j) = (Sw_w - Srw) * multx;
      aN (i, j) = (Sw_n - Srw) * multy;
      aS (i, j) = (Sw_s - Srw) * multy;
      aP (i, j) = aE (i, j) + aW (i, j) + aN (i, j) + aS (i, j);
    }
    

  /* ----- Boundary conditions ----- */

  applyBoundaryConditions2D();

    // I'm using an explicit scheme, so the boundary conditions are
    // different than in the implicit case.

  Range I(bi,ei), J(bj, ej);

    // Dirichlet right side
  aP(ei, J) += aE(ei, J);
  sp(ei, J) = 2 * aE(ei, J) * p(ei+1, J);
  aE(ei, J) = 0;

  // Neumann left side
  // Remember that Sw_e = 0.8 on the boundary, that is
  // the reason why sp(bi) has the form below.
  aP(bi, J) -= aW(bi, J) ;
  sp(bi, J) = aW(bi, J) * dx * ( 3.47e-7 * mu_w / k) ;
  aW(bi, J) = 0;
  
    return 0;
}

/*
 *  Lineal for realtive permeability , Upwind for Sw
 */
template<typename Tprec, int Dim>
inline bool BLES7<Tprec, Dim>::calcCoefficients3D () 
{
  static prec_t Sw_e, Sw_w, Sw_n, Sw_s, Sw_f, Sw_b;

    // Lineal 
  static prec_t multx = k * dt / 
    ( porosity * dx * dx * (1 - Srw - Sro) * mu_w );

  static prec_t multy = k * dt / 
    ( porosity * dy * dy * (1 - Srw - Sro) * mu_w );

  static prec_t multz = k * dt / 
    ( porosity * dz * dz * (1 - Srw - Sro) * mu_w );

  aE = 0.0; aW = 0.0; aN = 0.0; aS = 0.0; aF = 0.0; aB = 0.0; aP = 0.0; sp = 0.0;

  for (int ki = bk; ki <= ek; ++ki) 
    for (int i =  bi; i <= ei; ++i)
      for (int j = bj; j <= ej; ++j) {  

      // Upwind
	if ( phi_0(i+1, j, ki) >= phi_0(i, j, ki) ) Sw_e = phi_0(i+1, j, ki);
	else                                        Sw_e = phi_0(i  , j, ki);
	if ( phi_0(i-1, j, ki) >= phi_0(i, j, ki) ) Sw_w = phi_0(i-1, j, ki);
	else                                        Sw_w = phi_0(i  , j, ki);

	if ( phi_0(i, j+1, ki) >= phi_0(i, j, ki) ) Sw_n = phi_0(i, j+1, ki);
	else                                        Sw_n = phi_0(i, j, ki);
	if ( phi_0(i, j-1, ki) >= phi_0(i, j, ki) ) Sw_s = phi_0(i, j-1, ki);
	else                                        Sw_s = phi_0(i, j, ki);

	if ( phi_0(i, j, ki+1) >= phi_0(i, j, ki) ) Sw_f = phi_0(i, j, ki+1);
	else                                        Sw_f = phi_0(i, j, ki);
	if ( phi_0(i, j, ki-1) >= phi_0(i, j, ki) ) Sw_b = phi_0(i, j, ki-1);
	else                                        Sw_b = phi_0(i, j, ki);

      // Lineal
	aE (i, j, ki) = (Sw_e - Srw) * multx;
	aW (i, j, ki) = (Sw_w - Srw) * multx;
	aN (i, j, ki) = (Sw_n - Srw) * multy;
	aS (i, j, ki) = (Sw_s - Srw) * multy;
	aF (i, j, ki) = (Sw_f - Srw) * multz;
	aB (i, j, ki) = (Sw_b - Srw) * multz;
	aP (i, j, ki) = aE (i, j, ki) + aW (i, j, ki) + 
	  aN (i, j, ki) + aS (i, j, ki) + aF (i, j, ki) + aB (i, j, ki);
    }    

  /* ----- Boundary conditions ----- */

  applyBoundaryConditions3D();

    // I'm using an explicit scheme, so the boundary conditions are
    // different than in the implicit case.

  Range I(bi, ei), J(bj, ej), K(bk, ek);

    // Dirichlet right side
  aP(ei, J, K) += aE(ei, J, K);
  sp(ei, J, K) = 2 * aE(ei, J, K) * p(ei+1, J, K);
  aE(ei, J, K) = 0;

  // Neumann left side
  // Remember that Sw_e = 0.8 on the boundary, that is
  // the reason why sp(bi) has the form below.
  aP(bi, J, K) -= aW(bi, J, K) ;
  sp(bi, J, K) = aW(bi, J, K) * dx * ( 3.47e-7 * mu_w / k) ;
  aW(bi, J, K) = 0;

    return 0;
}


} // Tuna namespace


#endif //_BLES7_H_

















